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Night sky brightness at sites from DMSP-OLS satellite 
measurements 
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ABSTRACT 

We apply the sky brightness modelling technique introduced and developed by 
Roy Garstang to high-resolution DMSP-OLS satellite measurements of upward artifi- 
cial light flux and to GTOPO30 digital elevation data in order to predict the brightness 
distribution of the night sky at a given site in the primary astronomical photometric 
bands for a range of atmospheric aerosol contents. This method, based on global data 
and accounting for elevation, Earth curvature and mountain screening, allows the eval- 
uation of sky glow conditions over the entire sky for any site in the World, to evaluate 
its evolution, to disentangle the contribution of individual sources in the surrounding 
territory, and to identify main contributing sources. Sky brightness, naked eye stel- 
lar visibility and telescope limiting magnitude are produced as 3-dimensional arrays 
whose axes are the position on the sky and the atmospheric clarity. We compared our 
results to available measurements. 
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1 INTRODUCTION 

The change in the light in the night environment due to the 
introduction of artificial light is a true pollution, a grow- 
ing adverse impact on the night. Pollution means "impair- 
ment or alteration of the purity of the environment" or of 
its chemical/physical parameters. This alteration of natural 
light at night, called light pollution, can and does impact 
the environment and the health of the beings living in it 
(animals, plants and man), as shown by hundreds of sci- 
entific studies and reports (see e.g. Rich & Longcore 2002, 
Erren & Piekarski 2002, Cinzano 1994). The growth of the 
night sky brightness is one of the many effects of artificial 
light being wasted in the environment. It is a serious prob- 
lem. It endangers not only astronomical observations but 
also the perception of the Universe around us (see Craw- 
ford (1991), Kovalewski (1992), McNally (1994), Isobe & Hi- 
rayama (1998), Cinzano (2000a), Cohen & Sullivan (2001), 
Cinzano (2002), Schwarz (2003) and the International Dark- 
Sky Association Web site, www.darksky.org). The starry sky 
constitutes mankind's the only window to the universe be- 
yond the Earth. A fundamental heritage for the culture, 
both humanistic and scientific, and an important part of 
the our nighttime landscape patrimony is going to be lost, 
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for those alive today and for our children and their children. 
The worldwide growing interest for light pollution and its 
effects requires methods for monitoring this situation. 

The modelling of the brightness distribution of the night 
sky at a site is important to evaluate its suitability for astro- 
nomical observations, to quantify its sky glow, and to rec- 
ognize endangered parts of the sky hemisphere. Night sky 
models are useful in studying sky glow relationships with 
atmospheric conditions and to evaluate future changes in 
sky glow. The modelling is also required to disentangle the 
contribution of sources, such as individual cities, in order to 
recognize those areas producing the strongest impact and to 
undertake actions to limit light pollution. 

In 1986 Roy Garstang introduced a modelling tech- 
nique, developed and refined in the subsequent years 
(Garstang 1986, 1987, 1988, 1989a, b, 1991a, b, c, 1992, 
1993, 2000a), to compute the light pollution propagation 
in the atmosphere. He estimated the night sky brightness 
at many sites based on the geographical positions, altitudes 
and populations of the polluting cities. Cinzano (2000b) used 
Garstang models to disentangle the impact of individual 
cities constraining free functions with the condition that the 
sum of all the contributions with the natural sky brightness 
fits the observed sky brightness. However updated popu- 
lation data are not easily available worldwide, the upward 
light emission is not strictly proportional to the population. 
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Some polluting sources, such as industrial areas and air- 
ports, have very low population density but very high light 
emission. The U.S. Air Force Defence Meteorological Satel- 
lite Program (DMSP), Operational Linescan System (OLS) 
acquires direct observations of nocturnal lighting, making it 
possible to map the spatial distribution of nighttime lights 
(Sullivan 1989, 1991; Elvidge et al. 1997a, b, c, 2001, 2003a, 
b; Gallo et al. 2003; Henderson et al. 2003). Most nighttime 
OLS observations of urban centers are saturated, making 
the data of limited value for modeling purposes. However, 
Elvidge et al. (1999) were able to produce a radiance cali- 
brated global nighttime lights product using OLS data ac- 
quired at reduced gain settings, suitable for the quantita- 
tive measurement of the upward light emission (e.g. Isobe & 
Hamamura 2000, Luginbuhl 2001, Osman et al. 2001) and 
the evaluation of the artificial sky brightness produced by it 
(e.g. Falchi 1998; Falchi & Cinzano 2000). 

Cinzano et al. (2000) presented a method to map the 
artificial sky brightness across large territories in a given 
direction of the sky by evaluating the upward light emis- 
sion from DMSP high resolution radiance calibrated data 
(Elvidge et al. 1999) and the propagation of light pollu- 
tion with Garstang models. A World Atlas of the artificial 
night sky brightness at sea level was thus obtained (Cinzano, 
Falchi & Elvidge 2001b). This method was extended by Cin- 
zano, Falchi & Elvidge (2001a) to the mapping of naked 
eye and telescopic limiting magnitude based on the Schaefer 
(1990) and Garstang (2000b) approach and the GTOPO30 
elevation data. We extend and apply their method to the 
computation of the distribution of the night sky brightness 
and the limiting magnitude over the entire sky at any site 
for a range of atmospheric conditions and accounting for 
mountain screening. In sec. |2]we describe the computation 
on 3-dimensional arrays whose axes are the position on the 
sky and the atmospheric clarity and present our improve- 
ments. In sec. |H] we describe input data. In sec. 31 we deal 
with the disentangling of individual sources. In sec. [S] we 
discuss the application and in sec. |S|we present comparisons 
with available measurements. Conclusions are in sec. [7| 



2 COMPUTATION OF THE HYPERMAPS 

Artificial and natural sky brightness varies depending on the 
aerosol content of the atmosphere. The stellar extinction also 
vary substantially depending on the aerosol content of the 
local atmosphere. This in turn affects the limiting magni- 
tude. So any map of the sky of a site is a function of the 
aerosol content for which it has been computed. 

We refer to a hyper-map as a set of maps of the night 
sky brightness for a range of aerosol contents, b(z,ui,K), 
where z is the zenith distance, u> is the azimuth and K 
is the aerosol content expressed by the atmospheric clarity 
(Garstang 1986, 1989). As fig. Qshows, values on planes of 
the space of the variables perpendicular to the K axis give 
maps of the sky brightness for the given atmospheric clarity, 
values along a line parallel to the K axis give the brightness 
in the given point of the sky when the atmospheric aerosol 
content change, values along lines perpendicular to the K 
and u> axis give the sky brightness along an almucantar for 
the given atmospheric clarity. 

At a site in (x' , y') the hyper-map is given by: 



b(z,oj,K) 



e(x,y)f(x,y,x',y',z,u>,K) dx dy, (1) 



where f(x,y,x',y',z,u,K) is the light pollution propaga- 
tion function, i.e. the artificial sky brightness at (x' , y') in 
the direction given by (z,ui) per unit of upward light emis- 
sion e(x, y) produced by the unitary area in (x, y) when 
atmospheric aerosol content is K. If we divide a territory in 
land areas (h, I) with position (xh,yi), the hyper-map can 
be expressed as a tridimensional array bi t j t h given by: 



(2) 



bi,j,k - 22z-^ eh • l f( Xh ' yl,x ' ,y ' ,Zi,UJ: > , Kk ^ ' 

h I 

where en,i is the upward flux emitted by the land area 
(h,l), f(xh,yi,x',y',Zi,ujj,Kk) is the propagation function, 
Zi , ujj , Kk are an adequate discretization of the variables 
z,u,K and the sommatories are extended at all the land 
areas around the site inside a distance for which their con- 
tributions are non negligible. We divided the territory in 
the same land areas covered by pixels of the satellite data. 
We obtained the propagation function /, expressed as total 
flux per unit area of the telescope per unit solid angle per 
unit total upward light emission, with models for the light 
propagation in the atmosphere based on Garstang models 
(Garstang 1986, 1989): 

/=/(/3 m (/ l )/ m (w)+/?a(ft)/ a M)(H-D s )i(V,s)fi(u)du,(3) 

Jug 

where /3 m (/i) /3 a (/i) are respectively the scattering cross sec- 
tions of molecules and aerosols per unit volume at the alti- 
tude h, depending on the distance u along the line of sight 
of the observer, / m and / a are their normalized angular scat- 
tering functions (see sec !3.3ll . zu is the scattering angle, (,i(u) 
is the extinction of the light along its path from the scatter- 
ing volume to the telescope, i(ip, s) is the direct illuminance 
per unit flux produced by each source on each infinitesimal 
volume of atmosphere along the line-of-sight and (1 + Ds) 
is a correction factor which take into account the illumi- 
nance due at light already scattered once from molecules 
and aerosols which can be evaluated as Garstang (1984, 
1986), neglecting third and higher order scattering which 
can be significant for optical thickness higher than about 
0.5. Geometric relations and formulae accounting for Earth 
curvature have been given and discussed by Garstang (1989, 
sec. 2.2-2.5, eqs. 4-24). In Garstang's formulae the molecular 
scattering cross section per unit volume is /3 m = 2V m <7R. 

As Garstang (1989) and differently from Cinzano et al. 
(2001a) we take into account the elevation both of the source 
and of the site. 

Screening by terrain elevation was accounted as de- 
scribed in Cinzano et al. (2001a). The illuminance per unit 
flux was set in eq. (J^J to: 



»(iM=JM6/« a 



(4) 



where there is no screening by Earth curvature or by ter- 
rain elevation and i(lp,s) = elsewhere. Here I{ip) is the 
normalized emission function giving the relative light flux 
per unit solid angle emitted by each land area at the zenith 
distance ip, s is the distance between the source and the 
considered infinitesimal volume of atmosphere and £2 is the 
extinction along the light path, given by Garstang (1989). 
We check each point along the line of sight to determine if 
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the source area is blocked by terrain elevation or not, taking 
into account Earth curvature, by determining the position 
of the foot of the vertical of the considered point. Then we 
computed for every land area crossed by the line connecting 
this foot and the source area, the quantity cot ip (Cinzano 
et al. 2001a): 



cot ip — 



(A + E)- (h + E) cos(D/E) 
(h + E) sm(D/E) 



(5) 



where A is the elevation of the land area, D is the distance 
of its center from the center of the source area and E is the 
Earth radius. From it we determined the screening elevation 



h s = 



A + E 



cos(D*/E) ~ max(- cot ip) sin(D*/E) 



E. 



(6) 



where D* is the distance between the source area and the 
foot of the vertical, and h B is computed over the sea level. 
The illuminance i in eq. is set to zero when the elevation 
of the considered point is lower than the screening elevation. 
To speed up the calculation we computed only once the ar- 
ray, which gives the screening elevation for each point along 
the line of sight, for each azimuth of the line of sight and 
for each source, and we used it for any computation with 
different atmospheric parameters. We considered land areas 
as point sources located in their centres except when i = h, 
j — k in which case we used a four points approximation 
(Abramowitz & Stegun 1964). We assumed that the eleva- 
tion given by the GTOPO30 be the same everywhere inside 
each pixel. 

Another array was obtained for the natural sky bright- 
ness with the model introduced by Garstang (1989, sec. 3). 
The array &Ni,j,fc is the sum of (i) the directly transmitted 
light bd which arrives to the observer after extinction along 
the line of sight (Garstang 1989, eq. 30), (ii) the light scat- 
tered by molecules by Rayleigh scattering fo r (Garstang 1989, 
eq. 37), (iii) the light scattered by aerosols fe a (Garstang 
1989, eq. 32) : 



6n 



i,j,k 



= bdi,j,k + bri,j,k + b. 



'a,i,j,k ■ 



(7) 



In the computation of the natural sky brightness outside the 
scattering and absorbing layers of the atmosphere (Garstang 
1989, eq.29) we assumed as independent variables the bright- 
ness of a layer at infinity, due mainly to integrated star light, 
diffused galactic light and zodiacal light, and the brightness 
of the van Rhijn layer, due to airglow emission. 

The array of the total sky brightness is brij,k — hj,k + 
bNi,j,h' The sky brightness in the chosen photometric band 
was expressed as photon radiance (in ph cm -2 s _1 sr _1 ) or 
in magnitudes per arcsec 2 (Garstang 1989, eqs. 28, 39). 

We determined the observer's horizon computing the 
altitudes below which the line-of-sight encounter a screen 
by terrain, like e.g. a mountain, and set the total brightness 
to be zero below them. They are obtained evaluating the 
elevation /it of terrain at distance d along each azimuth 
direction and computing the maximum screening altitude 
angle i?: 



i? = maxarctan(/ix/rf) ■ 



(8) 



From the array of the total sky brightness in V band 
we can obtain a family of other arrays giving the naked-eye 
star visibility and the telescopic limiting magnitude. The 



magnitude over the atmosphere of a star at the threshold 
of visibility of an observer when the brightness of observed 
background is 6 b s in nanolambert and the stimulus size, 
i.e. the seeing disk diameter, is 6 in arcmin, has been given 
by Garstang (2000b) based on measurements of Blackwell 
(1946) and Knoll, Tousey & Hulburt(1946) and on a thresh- 
old criterion of 98 per cent probability of detection: 



"^star 

with: 

./ 

H = 

./ 

Ft ■■ 
F 2 -■ 



-13.98 -2.5 log «&/(<!+£). 



Fici(l + fci6 1/2 ) 2 (l + Q!<9 2 + yibHJ 2 ) , 
F 2 c 2 (l + k 2 b 1/2 ) 2 (l + a 2 9 2 + yib^O 2 ) ■ 
F a ,iFsc,iiwF c ,iF s ,i , 
F,, 2 Fsc,2F CB oF c oF B , 2 , 



(9) 

(10) 

(11) 
(12) 

(13) 



where i[ and i' 2 are the illuminations produced by the star, 
related respectively to the thresholds of scotopic and pho- 
topic vision, and the fraction is an artifact introduced by 
Garstang in order to put together smoothly the two compo- 
nents obtaining the best fit with cited measurements. Here 
F a takes into account the ratio between pupil areas of the 
observer and the pupil diameter used by the average of the 
Knoll, Tousey, Hulburt and Blackwell observers, Fsc takes 
into account the Stiles-Crawford effect, due to the decreas- 
ing of the efficiency to detect photons with the distance from 
the center of the pupil, producing a non linearity in the in- 
crease of sensibility when the eye pupil increase, F cs allows 
for the difference in color between the laboratory sources 
used in determining the relationships between i and b and 
the observed star, F D allows for star light extinction in the 
terrestrial atmosphere because star magnitudes are given 
outside the atmosphere, F s allows for the acuity of any par- 
ticular observer, defined so that F s < 1 implies an eye sen- 
sitivity higher than average due possibly to above average 
retinal sensitivity, scientific experience or an above average 
eye pupil size. Formulae have been given by Schaefer (1990) 
and Garstang (2000b) and applied by Cinzano et al. (2001a, 
eqs. 28-31) to which we refer the reader. The constants c, k, 
a, y, z in eq. 1101 are given by Garstang (2000b). The per- 
ceived background b ^ B is related to the total sky brightness 
under the atmosphere in V band given by our hyper-maps, 
converted from photon radiance to nanolambert (Garstang 
2000b): 



6 T /(F a F sc F cb ) . 



(14) 



where F c b allow for the difference in color between the lab- 
oratory sources and the night sky background, and F a , Fsc 
were already described. As a result we obtain the array nii t j t k 
of the visual limiting magnitude. The array of the telescopic 
limiting magnitudes can be calculated for the chosen instru- 
mental setup in a similar way (see the cited authors). 



3 INPUT DATA 

We summarize here the required input data, which has been 
already described and discussed by Cinzano et al. (2000, 
2001a). We refer the reader to their paper for details. We 
extended the input data to other continents in the same 
way. 
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3.1 Upward light emission data 

To compute the illuminance per unit flux i in eq.^Jwe need 
the relative intensity I(x, y, ij), \) emitted by every land area 
in (x, y) at azimuth \ arm zenith distance ip, i.e. the nor- 
malized emission function obtained measuring the relative 
emitted flux per unit solid angle per unit area in the di- 
rection ip and normalizing its integral to unity. If the land 
areas contain many light installations randomly distributed 
in type and orientation, we can assume this function axysim- 
metric I(x,y,ip). The corresponding absolute intensity is: 



I\x,y,i>) = e(x,y) x I(x,y,tp). 



(15) 



where e(x, y) is the total upward flux obtained from radiance 
calibrated data (Cinzano 2001a, eq. 35). 

We obtained the upward flux e(x, y) on a 30" x30" pixel 
size grid from the Operational Linescan System (OLS) car- 
ried by the DMSP satellites after a special requests to the 
U.S. Air Force made by the U.S. Department of Commerce, 
NOAA National Geophysical Data Center (NGDC), which 
serves as the archive for the DMSP and develops night time 
lights processing algorithms and products. OLS is an oscil- 
lating scan radiometer with low-light visible and thermal 
infrared (TIR) high-resolution imaging capabilities (Licskc 
1981). The OLS Photo Multiplier Tube (PMT) detector has 
a broad spectral response covering the range for primary 
emissions from the most widely used lamps for external 
lighting. The primary reduction steps were (Cinzano et al. 
2000, 2001a; Elvidge et al. 1999): 

(i) 1) acquisition of special OLS-PMT data at a number 
of reduced gain settings to avoid saturation on major ur- 
ban centres and, in the same time, overcome PMT dynamic 
range limitation. On board algorithms which adjust the vis- 
ible band gain were disabled. 

(ii) establishment of a reference grid with finer spatial 
resolution than the input imagery; 

(iii) identification of the cloud free section of each orbit 
based on OLS-TIR data; 

(iv) identification of lights, removal of noise and solar 
glare, cleaning of defective scan lines; 

(v) projection of the lights from cloud-free areas from 
each orbit into the reference grid; 

(vi) calibration to radiance units using preflight calibra- 
tion of digital numbers for given input telescope illuminance 
and gain settings in simulated space conditions; 

(vii) tallying of the total number of light detections in 
each grid cell and calculation of the average radiance value; 

(viii) filtering images based on frequency of detection to 
remove ephemeral events; 

(ix) transformation into latitude/longitude projection 
with 30"x30" pixel size; 

(x) Lucy-Richardson deconvolution to improve predic- 
tions for sites near sources (when possible this should be 
more properly done before step 7). 

(xi) Determination of the upward light intensity account- 
ing for the estimated atmospherical extinction in the light 
path from ground to the satellite, the assumed average spec- 
trum of night-time lighting (Cinzano et al. 2000a, eqs. 28-30) 
and the surface of the land areas. 

We can obtain I(x, y, ip) from the radiance measured in a 
set of individual orbit satellite images where the land area 



in (a;, y) is seen at different angles ip which are related to 
the distance D from the satellite nadir (Cinzano et al. 2000, 
eq. 17, 18). The emitted flux per solid angle per unit area 
in the direction ip is obtained from the measured radiance 
dividing by the extinction coefficient ^(ip) computed for 
curved-Earth (Cinzano et al. 2000, eq. 19). A study to ob- 
tain I(x, y, <p) in this way for every land area from DMSP- 
OLS individual orbit data is in progress (Cinzano, Falchi, 
Elvidge, in prep.). To be simple we assumed here that all 
land areas have in average the same normalized emission 
funct ion, g iven b y the parametric representation of Garstang 
Garstang (1986) in eq. (15) of Cinzano et al. (2000), which 
has been tested by studying in a single orbit satellite image 
the relation between the upward flux per unit solid angle 
per inhabitant of a large number of cities and their distance 
from the satellite nadir (Cinzano et al. 2000) and with many 
comparisons between model predictions and measurements 
by Garstang and by Cinzano (2000b). Likely it cannot be 
applied in areas where effective laws against light pollution 
are enforced or with unusual lighting habits. 

3.2 Elevation data 

As input elevation data we used GTOPO30, a global digi- 
tal elevation model by the U.S. Geological Survey's EROS 
Data Center (Gesch et al. 1999). This global data set covers 
the full extent of latitude and longitude with an horizontal 
grid spacing of 30" as our composite satellite image. The 
vertical units represent elevation in meters above mean sea 
level which ranges from -407 to 8,752 meters. We reassigned 
a value of zero to ocean areas, masked as "no data" with a 
value of -9999, and to altitudes under sea level. 



3.3 Atmospheric data 

In order to evaluate scattering and extinction we need a set 
of functions giving, for each triplet of longitude, latitude and 
elevation (x,y,h), the molecular and aerosol cross scatter- 
ing coefficients per units volume of atmosphere f3m(x,y,h) 
and /3 a (x,y,h), and the aerosol angular scattering func- 
tion f a (cu,x,y,h). The molecular angular scattering func- 
tion f m (u>) is known because it is Rayleigh scattering. The 
atmospheric data need to refer at a typical clean night in 
the chosen time of the year and to include information on 
denser aerosol layers, volcanic dust and Ozone layer. 

To be simple we applied here the standard atmospheric 
model already adopted by Garstang (1986, 1989) and Cin- 
zano et al. (2000, 2001a), neglecting geographical gradients 
and local particularities. It assumes: 

(i) the molecular atmospher e in hydro s tatic equilibrium 
under the gravitational force as lGarstand l|1986). 

(ii) the atmospheric haze aeros o ls nu mber density de- 
creasing exponentially as IGarstana <ll986T) 

(iii) a neglegible presence of sporadic denser aerosol lay- 
ers, volcanic dust and Ozone layer (studied by Garstang 
1991a, c). 

(iv) the normalized angula r scattering functi on for atmo- 
spheric haze aerosols given in IGarstana Ill991af) . 

(v) the aerosol content given by an atmospheric clarity 
parameter which measures the relative importance of aerosol 
and molecules for scattering light. 
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The Garstang atmospheric clarity parameter K measures 
the relative importance of aerosol and molecules for scatter- 
ing light in V band at ground level (Garstang 1996): 



K = 



&,. 



/3 m ,oll.lle- c » 



(16) 



where H is the altitude of the ground level over sea level 
and c is the inverse scale height of molecules. It assumes 
that there is only one ground level where all the polluting 
sources lie. To be more self-consistent when there are many 
cities at different elevations over sea level, we introduced an 
atmospheric clarity parameter K' defined at sea level: 



K' 



/3a,0 



/3m,0H.H 



(17) 



so that at ground level of each city K = K' &- c ~ a ' H , where a 
is the inverse scale height of aerosols. We can associate the 
atmospheric clarity K with vertical extinction (e.g. Garstang 
1991, eq. 6) and with other observable quantities like the 
horizontal visibility (Garstang 1989, eq. 38), the optical 
thickness r (Garstang 1986, eq. 22) and the Linke turbidity 
factor for total solar radiation (Garstang 1988). Extinction 
along light paths for this atmospheric model was given by 
Garstang (1989, eqs. 18-22). 



area in (xh,yi), searching for main polluting sources and 
making some statistic on their geographic distribution: 



bi,j,k{x h ,yi) = e hi if(x h ,yi,x',y',Zi,ujj,K k ) . 



(18) 



We can obtain hypermaps of sky brightness produced by 
each city or territory identifying pixels belonging to each city 
or territory of a given list and summing their contributions: 



h,j,k( n ) = / 

n-th city 



e-h,if{x h ,yi,x',y',z t ,u)j,K k ) . 



(19) 



Their comparison is helpful e.g. to understand if larger con- 
tributions come from few main cities or from many small 
towns, even in relation with atmospheric conditions. The 
fraction of sky brightness produced in a given array cell 
(i,j, k) by the sources inside a circular area of radius d can 
be obtained summing all contributions of land areas inside 
the distance d from the site and dividing by the sum of all 
contributions: 



l >*(d) = T > eh,if(.Xh,yi,x',y',Zi,Uj,K k ) ■ 

Vi.j.k * 'h,l 



(20) 



(x h -x') + 

(yi-y'?<id 2 



This is useful e.g. to evaluate the effectiveness of protection 
areas (Cinzano 2000c). 



3.4 Natural night sky brightness data 

The brightness bsi,j, due to integrated star light, diffused 
galactic light and zodiacal light, depends on the observed 
area of the sky and on the time. This dependence on the 
position of the sky is important when sky maps are made 
to quantify the visibility of astronomical phenomena, other- 
wise we can assume bsij constant and given by its average 
value in the considered site. The brightness of the Van Rhijn 
layer, &vr, depends on some factors like the geographical po- 
sition, the solar activity in the previous day, and the time 
after twilight. We referred our predictions to some hours 
after the twilight, when the night brightness decay at a con- 
stant value (Walker 1988, but see also Patat 2003a), and 
to minimum solar activity. If requested, the solar activity 
can be roughly accounted as Cinzano et al (2001a) or, more 
accurately, based on the correlation with the 10.7 cm solar 
radio flux (Walker 1988, Krisciunas 1999). The dependence 
of &vr by the geographical position suggests to study the 
natural sky brightness in the nearest unpolluted site, which 
can be located in the world atlas of artificial sky brightness 
(Cinzano et al. 2001b), in order to obtain bsi,j and 6vr- 
When only one or few measurements were available we as- 
sumed as Garstang (1989) bsi,j — OAbo and 6vr = 0.6&o 
and determined feo- 



4 DISENTANGLING INDIVIDUAL 
CONTRIBUTIONS 

We can make some analysis on the contributions from each 
30" x30" land area which enter in the summatory of eq.(2). 
First we can make hypermaps of sky brightness produced by 
individual land areas and compare them. Moreover, chosen 
an array cell of index (i, j, k) we can obtain a geographic map 
showing the contribution bij t k(xh, yi) produced by each land 



5 APPLICATION 

The software package lpskymap, written in Fortran-77, cal- 
culates the artificial night sky brightness, the total night 
sky brightness and the star visibility (limiting magnitude) 
over the entire sky at any site in the World. The availabil- 
ity of OLS-DMSP fixed gain data on a yearly or sub-yearly 
timescale will allow a fine time resolution. 

Results are arrays of the artificial night sky brightness, 
the total night sky brightness, the visual limiting magni- 
tude and the loss of visual limiting magnitude. Each hyper- 
map array is composed by a series of 19 x 37 pixel images in 
cartesian coordinates, one for each aerosol content K, spline 
interpolated over 91x181 pixels in cartesian coordinates or 
projected in 721x721 pixels in polar coordinates. Images go 
from to 360 degrees in azimuth, starting from East (in or- 
der to avoid to place the meridian at borders) toward South, 
and from horizon to zenith in altitude. They are saved in 
16-bit standard fits format with fitsio Fortran-77 routines 
developed by HEASARC at the NASA/GSFC. ASCII data 
tables are also provided. The night sky brightness in the cho- 
sen photometric band is given as photon radiance in ph s~ 2 
m~ 2 sr _1 or as astronomical brightness in mag arcsec~ 2 . 
Brightness in V band can be also expressed as luminance 
in /icd m~ 2 , using Garstang's conversion (Garstang 2002; 
Cinzano 2004). From the hyper- map arrays we can obtain: 

(i) sections perpendicular to the K axis: b(z,u),K — Kq). 
They are the maps of the sky brightness or limiting magni- 
tude for a given aerosol content and they correspond to each 
individual image of the series. 

(ii) secants parallel to the K axis: b(z — z ,u> — u>o,K). 
They provide the brightness or the limiting magnitude in a 
given point of the sky as the aerosol content changes. 

(iii) secants perpendicular to the K and u) axis: b(z, uj — 
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wo, K = Ko). They give the brightness or the limiting mag- 
nitude along an almucantar, e.g. the meridian, for a given 
aerosol content. 

The arrays computation steps are: 

(i) An input file is prepared with the geographical posi- 
tion and elevation of the site, the names of input DEM and 
lights frames and the position of their upper left corner. 

(ii) The array, i.e. the series of images, of the artificial 
night sky brightness is computed with the program lp- 
SKYMAP for a given range and step of the aerosol content, 
accounting for Earth curvature and elevation but not for 
screening. The radius of the contributing area can be 250 
km for sites in urbanized areas or 350 km for dark sites. 

(iii) Subimages with DEM and lights data have been 
cropped from the original large scale frames with the pro- 
gram makefrac. We use fits or RAW images 701x701 
px in size to limit the requirements of RAM memory 
during screening computation. They are checked for rela- 
tive mismatches which can be corrected with the program 

MAKESHIFT. 

(iv) The screening angles for each direction of observation 
and for each area inside a given radius from the site are 
computed with the program makescreen. We limited the 
radius to 200 km to avoid too long computation time. The 
program writes the screening data of each site in 106 files 
for a total size of 20GB uncompressed. It also calculates the 
horizon line as seen by the site. DEM pixels very near to the 
site are divided in 11x11 sub pixels evaluated separately. 

(v) An array containing the screened brightness is com- 
puted with the program lpskyscreen when there are rea- 
sons to believe that screening is not negligible. 

(vi) The images of the screened brightness array are sub- 
tracted from the correspondent images of the sky brightness 
array, after properly rescaling, in order to obtain the array 
of the night sky brightness corrected for mountain screening. 

(vii) The array is calibrated with the program lpskycal 
based on pre-flight calibration at 1996-1997, or on Cinzano 
et al. (2001b) calibration at 1998-1999 made with Earth- 
based measurements, or on observations taken at the same 
site. Measurements of Cinzano et al. (2001a) fitted predic- 
tions based on the pre-flight calibration with a ^0.35 mag 
arcsec -2 and a shift Am=-0.28 mag arcsec -2 , likely mainly 
due to the growth of light pollution in the period between 
the observations and the satellite data acquisitions. The pro- 
gram adds the natural sky brightness, producing a series of 
calibrated maps of the total night sky brightness, interpo- 
lated or not, and the limiting magnitude. It also adds the 
horizon line. It does not account for the refraction of light 
by the atmosphere which could increase the brightness near 
the horizon toward very far cities. 

(viii) Maps in polar coordinates are obtained with the 
program lpskypolar. East is up, North at right. 

(ix) Maps are analyzed with FTOOLS developed by 
HEASARC at the NASA/GSFC. 

(x) Comparison with observations is made with the pro- 
gram lpskycompare. Measurements should be "under the 
atmosphere". Statistical analysis is made with the software 
mathematica of Wolfram Research. 

A number of utility programs completes the package. The 
computation time depends on the geographical behavior of 



the site, like the quantity of dark pixels, the quantity of 
nonzero elevation pixels, etc. As an example, the compu- 
tation of one element of the array (i.e. a single map for 
a given atmospheric content) for Sunrise Rock on a work- 
station with Xeon processor running at 1700 MHz required 
about 35 hours for lpskymap, 10 hours for makescreen and 
6 hours for lpskyscreen. However, the computation with 
lpskymap for the site in Padua required 80 hours, even if 
restricted inside a radius of 250 km, whereas the same com- 
putation for Serra La Nave required 18 hours only. 



6 RESULTS 

In this section we present a sample of results which can be 
obtained with our method and some comparisons with avail- 
able measurements. Specific studies are reserved for forth- 
coming papers. 

NGDC request for low and medium gain DMSP-OLS 
data used in this work was granted from U.S. Air Force 
for the darkest nights of lunar cycles in March 1996 and 
January-February 1997. More recent data sets taken in the 
period 1999-2003 are already at our disposal, but they are 
still under reduction and, before we are able to use them, we 
need to solve a number of problems in the analysis process 
(Cinzano, Falchi & Elvidge, in prep.). Pre-flight calibration 
of upward flux refers to 1996-1997, to the average lighting 
spectra of Cinzano et al. (2000) and to an average verti- 
cal extinction in V band at imaging time assumed to be 
Am = 0.33 mag. All results are computed for minimum so- 
lar activity and refers to some hours after twilight. We tuned 
the parameter bo to fit the zenith natural sky brightness for 
clean atmosphere measured by Cinzano et al. (2001a) at 
Isola del Giglio, Italy, V — 21.74 ± 0.06 mag arcsec -2 in 
V band for average solar activity and 200 m altitude over 
sea level. It agrees well with the average natural night sky 
brightness of 21.6 mag arcsec -2 measured by Patat (2003a) 
at ESO-Paranal. In facts, the sky become darker going to 
lower elevation over sea level due to larger extinction, even if 
this phenomena is limited by the increase of the light scat- 
tered from aerosols and molecules along the line of sight 
(Garstang 1989). Patat (2003) reported a large contribu- 
tion from zodiacal light, about 0.18 mag arcsec -2 , which 
justifies the fact that he finds the sky slightly more lumi- 
nous than expected. The algorithm of Patat (2003b) applied 
to VLT images excludes almost completely the stellar com- 
ponent whereas Cinzano et al. (2001a) excluded only stars 
fainter than 18th mag, but the expected difference is only 
~0.03 mag arcsec -2 . The "visual" natural night sky bright- 
ness should be obtained from the measured one adding the 
average stellar background produced by stars with magni- 
tude Js7 missed by the instrument or the analysis (Cinzano 
& Falchi 2004). This contribute is about -0.26 mag arcsec -2 
when stars down to magnitude 24 are missed. In our bright- 
ness predictions we did not correct the natural night sky 
brightness to the visual value. 

Fig- |U shows the night sky brightness at Sunrise Rock, 
a site located in Mojave National Preserve, California, USA 
(long W 115° 33'6.4",latN35° 18' 57.7") at 1534 m over sea 
level. This site is mainly polluted by the lights of Las Vegas, 
about 100 km away. Azimuth goes from to 360 degrees, 
starting from East toward South. Fig. [3] shows the night 
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sky brightness screened by mountains, which amounts to 
few hundredth of magnitude. Fig.^Jshows a comparison be- 
tween predictions for atmospheric clarities K'=0.5 (squares) 
or K'—3 (crosses) and Uband measurements taken on 2003 
May 8 at 05.34-06.00 UT (Duriscoe et al. 2004) with vertical 
extinction kv=0.18 mag. The agreement is excellent after an 
uniform scaling of about -0.3 mag arcsec -2 . It suggests an 
increase of light pollution from 1997 to 2003 of «5 per cent 
per year, slightly less than the average yearly growth ~6 
per cent estimated by Cinzano (2003). A comparison with 
a data set taken on 2003 September 22 at 06.27-06.58 UT 
with kv=0.26 mag gives similar results. 

Fig. shows the night sky brightness at Serra la Nave 
Observatory (long E 14° 58' 24", lat N 37° 41' 30") at 
1734 m over sea level on the Mt. Etna volcano, Italy. This 
site is situated at few kilometers from a densely populated 
area with ~1.8 10 6 inhabitants, which includes the cities of 
Catania (23 km) and Messina (75 km). Fig.[£]shows a com- 
parison between predictions for atmospheric clarities K'—l 
(squares) and K'=2 (crosses) with V band measurements 
taken on 1998 February 22-23 at 18.00-20.00 UT with ver- 
tical extinction kv=0.26 mag (Catanzaro & Catalano 2000; 
see fig. 2). The agreement is good. The fit is slightly bet- 
ter for the model with K'=l, corresponding to a vertical 
extinction of ky=0.17 mag, which is smaller than the mea- 
sured one. However the vertical extinction at this site could 
be locally determined by the volcanic dust (Catanzaro, priv. 
comm.) whereas K' depends on the average aerosol content 
of the entire area with 250 km radius, so they do not need 
to match. 

The effect of an increase of the aerosol content depends 
on the distribution of sources around the site. In general 
it decreases the zenith brightness when the distance of the 
main sources is larger than few kilometers, decreases the 
brightness at low elevation in direction of far sources, in- 
creases the brightness at very low elevation in direction of 
sources at small or average distance. This could explain the 
different behavior of the sky brightness with the aerosol con- 
tent in these two sites. 

Fig. Q shows the night sky brightness versus the zenith 
distance at G. Ruggeri Observatory, Padova, Italy (long E 
11° 53' 20", lat N 45° 25' 10"). This site is located inside 
a city of 8 10 5 inhabitants in a plain with more than 4 
10 6 inhabitants. Positive zenith distances collect measure- 
ments with azimuth inside ±90° from the direction of the 
city centre. Open symbols are the V band measurements 
taken on 1998 March 26 at 20.00-23.30 UT, with k v =0A8 
mag (Favero et al 2000). Filled symbols are predictions in the 
same directions for atmospheric clarity K'=3, correspond- 
ing to kv=0.65 mag. For smaller values of K' the brightness 
is underestimated by a constant value. This is likely due 
to the fact that our model cannot accurately account for 
the scattered light coming from lighting installations inside 
few hundreds of meters from the site because pixel sizes are 
of the order of 1 km. We used for this prediction the cali- 
bration made for 1998-1999 by Cinzano et al. (2001a). For 
an atmospheric clarity K ^ 2.2, i.e. for an optical depth 
t ^ 0.5, the double scattering approximation could be not 
fully adequate (Garstang 1989; Cinzano et al. 2000). Fig.QO 
shows the contribution to the artificial night sky brightness 
produced in the same site from the sources outside Padua 
for atmospheric clarity K =1.9 (ky=0.48 mag). The area 



neglected in the prediction is shown in Fig. [5] together with 
the distribution of lights in the Padana Plain surrounding 
Padua from OLS-DMSP satellite data. 

Fig. I1UI shows in polar coordinates the total 
night sky brightness in V band at Mt. Graham Ob- 
servatory, USA (long W 109° 53' 31", lat N 32° 
42' 5", 3191m o.s.l.) for atmospheric clarity K' = 
0.5. It can be compared with the image available 
at the web address http://mgpc3.as.arizona.edu/images/ 
Night%20Sky%201arge.jpg or with fig. 8 of Garstang (1989), 
which shows only the artificial brightness. 

Fig. Illl shows the naked eye limiting magnitude at Sun- 
rise Rock. Limiting magnitude is computed for observers of 
average experience and capability F s — 1, aged 40 years, 
98 per cent detection probability (faintest star that the ob- 
server sees surely and not the faintest suspected star) and 
star color index B — V — 0.7 mag. Experienced amateur as- 
tronomers are more trained in naked-eye observation than 
unexperienced peoples and can consider detected a star at 
a much smaller detection probability so that their limiting 
magnitude can be more than one magnitude larger (Schaefer 
1990). See the discussion in Cinzano et al. (2001a). 

We checked the effects of the mountain screening trying 
to reproduce the umbrae on the sky modelled by Schaefer 
(1988) and due to the screening produced by the Mauna 
Kea on the light of the rising sun backscattered to the ob- 
server. Fig. ll2l shows the analogous of the Schaefer's umbrae 
produced by a source of light pollution instead of the Sun. 
A city screened by a large conic mountain (left) projects 
an umbra over the horizon (right). When the mountain is 
off-set in respect to the line observer-source, a non sym- 
metric penumbra appears. Here the penumbra is at higher 
altitudes than in the Wynn- Williams' photo (Schaefer 1988, 
fig.l) likely because the observer is at lower elevation. Fur- 
ther examples of umbrae and baffles are shown by Cinzano 
& Elvidge (2003a fig.l and fig.3, 2003b). 



7 CONCLUSIONS 

We extended the seminal works of Garstang by applying his 
models to upward flux data from DMSP satellites and to 
GTOPO30 digital elevation models, and by accounting for 
mountain screening. The presented method allows one to 
monitor the artificial sky brightness and visual or telescopic 
limiting magnitudes at astronomical sites or in any other 
site in the World. 

This study provides fundamental information in evalu- 
ating observing sites suitable for astronomical observations, 
to quantify sky glow, to recognize endangered parts of the 
sky hemisphere when measurements are not readily avail- 
able or easy feasible, and to quantify the ability of the res- 
ident population to perceive the Universe they live in. The 
method enables to study the relationship of night sky bright- 
ness with aerosol content and to evaluate its changes with 
time. The method also allows one to analyze the adverse 
impacts on a site from the surrounding territories, making 
it possible to disentangle individual contributions in order 
to recognize those that are producing the stronger impact 
and hence to undertake actions to limit light pollution (the 
use of fully shielded fixtures, limitation of the downward 
flux wasted outside the lighted surface, use of lamps with 
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reduced scotopic emission, flux reduction whenever possi- 
ble, no lighting where not necessary, restraining of lighting 
growth rates or lighting density, etc). We also present some 
tests of the method. The effects of light pollution on the 
night sky are easily evident in the maps in the text. 

Important refinements needs to be done in the future 
years: i) it may be possible to derive the angular distribution 
of light emissions from major sources of nighttime lighting 
from OLS or future satellite data (Cinzano, Falchi, Elvidge 
in prep.). This will improve accuracy of the modelling, in 
particular where laws against light pollution are enforced; 
ii) a global Atlas of the growth rates of light pollution and 
zenith night sky brightness from satellite data (Cinzano, 
Falchi, Elvidge in prep.) will make it possible to predict the 
evolution of the night sky situation at sites; iii) a worldwide 
atmospheric data set giving the atmospheric conditions in 
any land area for the same nights of satellite measurements 
or for a typical local clear night will allow to replace the 
standard atmosphere with the true atmosphere or the typ- 
ical local atmosphere; iv) the availability of spectra of the 
light emission of each land area from satellite will allow a 
more precise conversion of OLS data to astronomical pho- 
tometrical bands and an accurate modelling of the colors of 
the night sky; v) a large number of accurate measurements 
of night sky brightness and visual limiting magnitude in- 
cluding the evaluation of the atmospheric content, from e.g. 
the vertical extinction, will allow to better constrain predic- 
tions allowing improvements of the modelling technique. The 
International Dark-Sky Association, the organization which 
takes care of the battle against light pollution and the pro- 
tection of the night sky is making a large worldwide effort to 
collect accurate measurements of both night sky brightness 
and stellar extinction (e.g. Cinzano & Falchi 2004). They 
constitute a fundamental component of the monitoring of 
the night sky situation in the World. 
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Figure 1. Projections of the hyper-map on different planes. 




Figure 2. Night sky brightness at Sunrise Rock, USA for atmospheric clarity K'=0.5. 




Figure 3. Brightness screened by mountains at Sunrise Rock, USA for atmospheric clarity K'=0.5. Each level from blue to violet is 
0.01 mag arcsec - 2 . 
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Figure 4. Comparison between predictions and V band measurements at Sunrise Rock for atmospheric clarities A"'=0.5 (squares) and 
K'=3 (crosses). Units are mag arcsec -2 . 
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Figure 5. Night sky brightness at Serra la Nave Observatory, Italy for atmospheric clarity K'=l. 
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Figure 6. Comparison between predictions and V band measurements at Serra la Nave Observatory for atmospheric clarities K'=l 
(squares) and K'=2 (crosses). Units are mag arcsec -2 . 
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Figure 7. Brightness - zenith distance relation measured at G. Ruggeri Observatory, Italy (open symbols) and predictions for the same 
viewing directions (filled symbols) for atmospheric clarity K'=3 versus the zenith distance. Positive elevations collect measurements with 
zenith distances less than ±90° from the direction of the city centre. 
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Figure 8. Contribution to the artificial night sky brightness at Padua from the sources outside Padua for atmospheric clarity K'=l. 




Figure 9. Distribution of lights in the plain surrounding Padua from OLS-DMSP satellite data. Dark section is the neglected area in 
the prediction of Fig. [S] The region shown is 50' square in geographic latitude/longitude projection (approximately 65 X 93 km). 
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Figure 10. Night sky brightness in V band at Mt. Graham Observatory, USA in polar coordinates for atmospheric clarity K' = 0.5. 
The figure is plotted with East at bottom, North at left. 
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Figure 11. Naked eye limiting magnitude at Sunrise Rock, USA for atmospheric clarity K' = 0.5 and 98 per cent detection probability. 




Figure 12. A city screened by a large mountain (left), off-set in respect to the line observer-source, projects an asymmetric Schaefer's 
umbra on the sky (right). Brightness scale is arbitrary. 
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